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ABSTRACT 

High-energy, multi-component plasmas in which pair creation and annihilation, lepton- 
lepton scattering, lepton-proton scattering, and Comptonization all contribute to establishing 
the particle and photon distributions, are present In a broad range of compact astrophysical 
objects. The different constituents are often not in equilibrium with each other, and this mixture 
of interacting particles and radiation can produce substantial deviations from a Maxwellian 
profile for the lepton distributions. Earlier work has included much of the microphysics needed to 
account for electron-photon and electron-proton interactions, but little has been done to handle 
the redistribution of the particles as a result of their Coulomb interaction with themselves. The 
most detailed analysis thus far for finding the exact electron distribution appears to have 
been done within the framework of non-thermal models, where the electron distribution is 
approximated as a thermal one at low energy with a non-thermal tall at higher energy. Recent 
attention, however, has been focused on thermal models. 

Our goal here is to use a Fokker-Planck approach in order to develop a fully self-consistent 
theory for the interaction of arbitrarily distributed particles and radiation to arrive at an 
accurate representation of the high-energy plasma in these sources. We conduct several tests 
representative of two dominant segments of parameter space. For high source compactness 
of the total radiation field, I ~ 10^, we find that although the electron distribution deviates 
substantially from a Maxwellian, the resulting photon spectra are insensitive to the exact 
electron shape, in accordance with some earlier results. For low source compactness, I ~ few, 
and an optical depth < 0.2, however, we find that both the electron distribution and the photon 
spectra differ strongly from what they would be In the case of a Maxwellian distribution. In 
addition, for all values of compactness, we find that different electron distributions lead to 
different positron number densities and proton equilibrium temperatures. This means that the 
ratio of radiation pressure to proton pressure is strongly dependent on the lepton distribution, 
which might lead to different configurations of hydrostatic equilibrium. This, in turn, may 
change the compactness, optical depth, and heating and cooling rates, and therefore lead to 
an additional change In the spectrum. An important result of our analysis, is the derivation 
of useful, approximate analytical forms for the electron distribution in the case of strongly 
non-Maxwelllan plasmas. 

Subject headings: accretion disks — black hole physics — galaxies: Seyfert — gamma rays: theory — 
plasmas — radiative transfer 
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1. Introduction 

There is a great deal of interest in the physics of relativistic, multi-component, astrophysical 
plasmas, which appear to be present in a wide variety of high-energy emitting objects, such as 
neutron stars, black holes, and Active Galactic Nuclei (AGNs). The particle distributions in these 
sources depend critically on key physical processes, such as pair creation and annihilation, lepton- 
lepton scattering, lepton-proton scattering, and Comptonization. Often, the particle population 
and the radiation field are not in equilibrium, and it is possible that in some cases, such as the 
inner coronal regions in black hole systems like Cygnus X-1 and 1E1740. 7-2942, the electrons and 
protons are themselves not equilibrated, with the proton "temperature" greatly exceeding that of 
the electrons and positrons (Shapiro, Lightman & Eardley 1976; for a more recent discussion and 
more extensive list of references, see also Misra & Melia 1996). 

The attention paid to these plasmas is fostered by continued, exciting observations of the 
underlying objects by a number of high-energy instruments on satellites such as the Compton 
Gamma-Ray Observatory (CGRO). For example, OSSE and COMPTEL on CGRO have recently 
revealed that persistent MeV gamma-rays are emitted by several Galactic black hole candidates 
(GBHCs), including Cygnus X-1 (Johnson et al. 1993; McConncU ct al. 1994) and GRO J0422+32 
(van Dijk et al. 1995). Single component thermal models (e.g., Sunyaev & Titarchuk 1980; 
Titarchuk 1994) have difficulty accounting for these high-energy tails, which strongly hint at 
nonthermal processes. In these systems, not only is a substantial fraction of the radiation produced 
externally and therefore not initially in thermal equilibrium with the electron-positron pair plasma, 
but the combination of pair creation/annihilation (Liang 1979; Kusunose 1987; Svensson 1990; 
Melia & Misra 1993; Misra & Melia 1996), thermal and non-thermal radiation, and the heating 
of protons via gravitational energy dissipation, results in a mixture of interacting particles and 
radiation that produces important deviations from a Maxwellian profile for the lepton distribution. 
Recently, Li, Kusunose Sz Liang (1996) discussed a stochastic particle acceleration scheme to 
account for this high-energy emission, and concluded that the energized particles consist of a 
blend of thermal and nonthermal components. In addition, they found that for certain parameter 
regimes, acceleration is more efficient than the cooling processes and must therefore be taken into 
account when determining the particle distributions of the emitting plasma. Most importantly, 
they showed that a deviation of the particle distribution from a Maxwellian at high energies can 
account for most of the emission above 1 MeV via inverse-Compton scattering. 

This alone makes the fully-consistent approach we describe here essential for a more complete 
understanding of the high-energy emission in GBHCs, but there are several equally important 
additional aspects to this problem. One of these is the appearance of a transient ~ McV bump 
discovered in Cygnus X-1 by Ling et al. (1987), which is difficult to explain with the preliminary 
study of Li et al. (1996). Earlier, Melia & Misra (1993) suggested that the transient bump might 
be due to bremsstrahlung-self-Compton processes within the gravitationally-compressed, inner 
hot region of the disk. However, without a complete analysis of the actual particle distribution, 
including its deviations at high energy, this hypothesis remains to be tested. A second important 
concern is the rapid spectral variability in these sources, which the Rossi X-Ray Timing Explorer is 
now beginning to sample with sub-millisecond resolution. Since the rapid acceleration and cooling 
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processes affect the particle energy re-distribution on this time scale (e.g., Melia & Misra 1993), 
it would be highly desirable to have a time-dependent scheme to track the particle and spectral 
evolution in a self-consistent manner. The approach described in this paper will be able to address 
each of these issues. 

The OSSE instrument on CGRO has also detected several Seyfert galaxies above 60 keV (e.g., 
Maisack ct al. 1993; Cameron et al. 1995; Madejski 1995), showing a break or an exponential 
cut-off at ~ 50 — 100 keV, with an overall spectrum not unlike those of the GBHCs (e.g., Sunyaev 
et al. 1991). Together with X-ray observations by GINGA, these OSSE data suggest that a single 
robust mechanism may be operating in both classes of sources (Fabian 1994). However, although 
basic models such as the pure non-thermal pair model (for a review, see Svensson 1994) predict 
a steepening of the hard X-ray spectrum, they do not easily produce a break as sharp as that 
seen by OSSE. In addition, they predict a flattening of the spectrum at ~ 100 — 200 kcV, and an 
annihilation line, neither of which has been observed so far. The Seyfert spectra may instead be 
based on thermal models, or more realistically, quasi-thermal models (e.g., Haardt Sz Maraschi 1991, 
1993; Ghisellini, Haardt & Fabian 1993; Titarchuck & Mastichiadis 1994; Zdziarski et al. 1994; 
Zdziarski et al. 1995; Haardt, Maraschi Sz Ghisellini 1996; see Svensson 1996 for a recent review). 
In the past, purely thermal models for the high-energy spectra of AGNs have been criticized on the 
basis of the long thermalization time needed to attain thermal equilibrium at high temperatures. 
The reason for this is that the observed X-ray variability time scale is shorter than the inferred time 
required for two-body thermalization, suggesting that the plasma cannot be Maxwellian. Thus, 
both thermal and non-thermal models for the X-ray spectra in Seyfert galaxies have significant 
problems. In some regions of parameter space, non-thermal and thermal particle distributions 
produce very similar spectra, as long as the radiative process is multi-Compton scattering and 
the maximum particle energy is a few MeV (Ghisellini et al. 1993), suggesting that a compromise 
model may still work. However, this similarity in spectra occurs only for a restricted range of optical 
depths and source photon compactness, and the assumed dynamics of the emitting plasma may not 
be consistent with the luminosities implied by the observed photon emissivities. There is clearly 
a need to establish the correct (non-Maxwellian) particle distribution in these sources for a broad 
range of system parameters. In a portion of the available phase space, the photons are not multiply- 
scattered so the resultant spectrum will in fact reflect the underlying particle distribution. Very 
importantly, the self-consistently determined particle and photon populations must be reconciled 
with the permitted system-dependent considerations, such as the allowed compactness and lepton 
temperature, which turn out to be dependent upon the actual electron and positron distributions 
for part of the parameter space. 

Our goal in this paper is to begin the process of developing a fully self-consistent theory for the 
interaction of arbitrarily distributed particles and radiation to arrive at an accurate representation 
of the high-energy plasma in compact sources. We first consider the non-Maxwellian parameter 
space in §2. We then motivate the use of the Fokker-Planck approach in §3, and we derive the 
general Fokker-Planck coefficients for arbitrary particle distributions in §4. In §5 we extend the 
generality of our scheme by including the electron/positron-proton interaction and the effects due to 
Comptonization (pair creation and annihilation are included by means of some standard formulae 
available in the literature and discussed in §6). We test the algorithm in §6. Much attention in that 
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section is paid to the physical consequences of the distribution function being non-MaxweUian. We 
end with our conclusions in §7. The Appendix describes in detail the numerical time-dependent 
scheme that has been used to solve the full electron-photon-proton set of equations. The reasons 
for neglecting large angle scatterings in the diffusion coefficient are also given in the Appendix. 

2. The Non-M£Lxwellian Parameter Space 

In this section, we discuss briefly the circumstances under which one should reasonably expect 
to see the electron distribution deviating substantially from a Maxwellian shape. Although this 

issue is covered extensively in the literature (e.g., Stepney 1983; Guilbert and Stepney 1985; DL89; 
Baring 1987; Ghisellini, Haardt & Fabian 1993, hereafter GHF; Fabian 1994), there seems to be no 
simple expression which would characterize the parameter space where the non-thermal distribution 
can form for a general source geometry and arbitrary optical depth. Our goal here is to provide 
such an expression. 

The most likely process to dominate the truncation of the electron distribution is inverse 
Compton cooling on low energy photons (Fabian 1994). As is well known (e.g., Rybicki & Lightman 
1979), the rate of energy transfer in the Thomson scattering regime, from a single electron with 
gamma-factor 7 to soft photons with energy density Uph, is 

mec^ (^)" = ^ arc ^VC/ph, (1) 

where Up^ is the soft photon energy density within the source. For a source with a spherical 
geometry of radius R, one can approximate this expression as 

where / = Lgt / ^T^Rm^c^ is the compactness of the region, and L is the total source luminosity. 
Equation (2) does not change for a disk geometry, but the compactness I should then be replaced 
by I = Lar / hrrieC^ , where h is the geometrical thickness of the disk, and L is the luminosity of a 
cube (i.e., the standard deflnition of I; e.g., Svensson 1996; § 6). 

For the electron-electron energy exchange rate, we will use an approximation that is well known 
in the literature and re-derived below (Eq. 30). This approximation is valid for E (s> &e) equal 
to the average thermal 7-factor (7) : 

{d_E\ ^_(3/2)i^iA_i_. (3) 



\ dt / ee tT 

Let us now introduce the parameter A which will give the ratio of the inverse Compton cooling 
rate to the e-e Coulomb energy transfer rate for the average electron: 



A = 



;f )..l ^ '-rlnA 

where tt = n^cTTR is the Thomson optical depth. Considering the relativistic and non-relativistic 
limits of (Pj)^, and merging these two expressions for the mildly relativistic domain, one can 
simpliiy Equation (4) to 
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The physical significance of this parameter is such that when A is of order unity, one should 
expect a noticeable truncation in the tail of the Maxwellian distribution. It is not a quantitative 
parameter that describes exactly how the high energy tail of the Maxwellian is truncated, because 
this is determined not only by electron-photon and electron-electron interactions, but also by the 
interaction that heats the electrons. Note that A ^ 1 does not necessarily mean that the photon 
spectrum will be changed significantly due to such a truncation (see the tests in § 6 below) . 

One can rewrite Equation (5) in terms of the maximum compactness of the region where e-e 
Coulomb collisions can still maintain a Maxwellian distribution. This gives an expression similar 
to the one derived by Fabian (1994), except for the factor (1 -|- tt)/tt, which was approximated 
by unity for tt ~ 1. Since recently much attention has been given to optically thin plasmas (e.g., 
Svensson 1996), we will here retain this factor. 

Another competitive thermalization process is the so called synchrotron boiler (e.g., Ghisellini, 
Guilbert &; Svensson 1988). It turns out that synchrotron photons emitted and re-absorbed by 
electrons can be a very important mechanism for the electron thermalization. Although we do 
not consider this thermalization mechanism in this paper, i.e., we assume that magnetic field is 
small, it is necessary to asses the importance of this process. The total synchrotron power emitted 
by an electron is given by an equation similar to Equation (1), except that the radiation energy 
density should be replaced by the magnetic field energy density (J7mag = -B^/Btt). Thus, when the 
magnetic field energy density is larger than the radiation energy density, the electron distribution 
will be thermalized due to the synchrotron self-absorption process or due to Coulomb electron- 
electron interactions, if the latter dominate. To see when synchrotron thermalization dominates 
over Coulomb thermalization, let us define a 'magnetic' thermalization parameter, A^. Following 
Ghisellini, Guilbert & Svensson (1988), we first define the magnetic compactness parameter 



The synchrotron self-absorption will dominate as a thermafizing mechanism if A^, defined as the 
ratio of the synchrotron emissivity to the electron-electron energy exchange rate, is larger than 
unity: 



Here we have followed the same steps as those leading to Equation (5). The physical significance 
of Xm is that if it larger than unity, then synchrotron self-absorption thermalization dominates 
over Coulomb thermalization. It is interesting to note that recent models for the origin of X- 
rays in Seyfert Galaxies invoke short-lived and intense magnetic fiares on the surface of a cold 
accretion disk (e.g., Haardt, Ghisellini & Maraschi 1994; Nayakshin & Melia 1997a). In these 
fiares lb > I, since the conditions for plasma confinement require that radiation pressure (the 
dominant pressure in the case Z S> 1) is much smaller than the magnetic stress (Nayakshin &; Melia 
1997b). Therefore, it is very likely that the synchrotron processes keeps the particles thermal in 
such fiares, unless the heating mechanism acts in the opposite direction (the heating mechanism is 
likely to be of a collective rather than a two-body nature, and it is still unknown; see § 6 below). 
In other circumstances, e.g., inside of the accretion disk, the magnetic field energy density must 



(6) 



mag 




(7) 
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certainly be below the equipartition value, and for the radiation dominated part of the disk, the 
'synchrotron thermalization' will be weaker than inverse Compton cooling. For our purposes here, 
we restrict the parameter space of interest to the regions where is smaller than 1, so that 
Coulomb thermalization dominates. 

We should also note that other processes are unlikely to be of a real importance in determining 
the shape of the equilibrium electron distribution (this will be clear from the discussion below). 
For example, pair creation and annihilation is a considerably slower process; we experimented with 
non-thermal models where pairs are injected with some power-law distribution (e.g., Lightman 
k. Zdziarski 1987), and found that pairs have essentially the same distribution as the electrons 
everywhere, except that normalization of the power-law component is different with respect to the 
thermalized part of the distribution. Thus, it is fair to assume that the electrons and positrons have 
the same distributions but with different normalizations. Next, proton heating is also unlikely to 
change the shape of the electron distribution because it is strongest at the low energy end. In that 
region, the equilibrium distribution when only the proton-electron interactions remain in the full 
equation is a Maxwellian with the proton temperature. But for the low energy end {E <C kT), the 
electron distribution function is the same for all temperatures, differing only by its normalization. 
Therefore, the low energy end of the electron distribution function in the presence of proton heating 
has the Maxwellian shape, irrespective of how strong the proton heating is. If the electron heating 
mechanism is different from proton heating, however, then there is a possibility that it will lead to 
deviations from a Maxwellian on the low energy end (where it is the strongest, presumably), but 
we leave this for future work since the exact heating mechanism is not yet known. 

3. The Fokker-Planck Approximation For Coulomb Interactions 

As is well known, the Coulomb cross-section diverges for small angle scattering, where the 
relative change in the particle's energy approaches a negligible value (e.g., Chandrasekhar 1942). 
It is therefore unwise and numerically very challenging to treat Coulomb scattering using the 
otherwise commonly employed Boltzmann collision integral. Instead, this feature of the Coulomb 
force lends itself more naturally to a Fokker-Planck (FP) approach for solving the kinetic equation 
in the presence of an type of interaction. The FP procedure is particularly powerful when 
the change in particle energy /S.E is much smaller than its incident energy E. When this holds, 
the distribution functions appearing inside the Boltzmann collision integral can be decomposed as 
power series in the small expansion parameter /S.E/E (^ 1). An example of how this works in 
practice is the Kompaneets equation (Kompaneets 1957), which is well-known in astrophysics. 

When the particle distribution function is isotropic and homogeneous in space, the FP equation 
for the distribution function takes the form 



df{E,t) _ d 
dt dE L 



a{E,t)f{E,t) D{E,t)f{E,t) , (8) 



where f{E, t) is the distribution function of electrons in energy space (i.e., f{E, t) = Pj'^f*{E, t), E 
is kinetic energy of the electron in units of meC^, E = (j — 1), /3 = {1 — 1/7^)"^/^ and f{E,t)dE = 
f*{'p,t)d?p gives the volume density of electrons in the phase space element dPp, and p is the 
electron momentum in units nieC^), a{E,t) is the average energy exchange rate and D(E,t) is the 
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energy dispersion rate of the test particle. Both a{E, t) and D{E, t) are functions of time because 
the distribution function over which they are convolved depends on time. Since E is dimensionless, 
we will use E and (7— 1) interchangeably throughout the paper, and 7(7) is the same distribution 

as f{E). 

In their application of this equation to the study of electron thermalization in a background 
relativistic Maxwell-Boltzmann (MB) plasma, Dermer and Liang (1989) derived the FP coefficients 
for a small population of particles scattering oflF a relativistic MB distribution. In this limit 
where the test particle number density is much lower than that of the background MB plasma, 
Equation (8) is linear in the distribution function f{E,t). However, this approach fails when 
the real electron distribution differs substantially from the perfect (background) MB profile. The 
original motivation for taking this simplified approach was to produce onc-dimensional integrals 
for the FP coefficients. However, as we shall see in §4 below, the FP coefficients derived there 
for an arbitrary particle distribution are themselves simply one-dimensional integrals and there is 
therefore no computational advantage to be had in restricting the problem to the quasi-thermalized 
initial condition if the temperature is allowed to evolve. Of course, our generalized approach is 
practical only if we can find a numerical scheme that makes the problem tractable. The nature of 
the computational scheme that solves this problem is discussed in detail in the Appendix. 
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4. Coulomb Electron-Electron Fokker-Planck Coefficients 
For Arbitrary Particle Distributions 

4.1. The Single Particle Energy Exchange 0(7,71) and Diffusion D(7,7i) Coefficients 

We begin by deriving the Fokker-Planck coefficients for e-e Coulomb interactions for a mono- 
energetic particle distribution. The FP coefficients for an arbitrary particle profile can be calculated 
from these by convolving them over f{E,t). We consider the FP energy exchange and diffusion 
coefficients for a test particle with dimensionless energy E = (7 — 1) (all electron energies are in 
units of TUeC^ in this paper) Coulomb scattering off a mono-energetic electron (positron) distribution 
with energy Ei = (71 — 1). Our notation and some of our starting equations are based on earlier 
work by Dermer (1985) and Dermer 8z Liang (1989) (DL89 hereafter). 

The distribution function is normalized by the condition 

1 = J djPj^r{j) = I d7/(7)- (9) 
1 1 

For two such functions /i (71) and /2 (72), the relativistically correct reaction rate is 

00 00 1 

r = 2(1+7x2) / / dlc{ll-l?" j dtxA*(7i)r(72)c7(7.,7c), (10) 

11 -1 

where ni and n2 are the corresponding particle number densities, 7^. is the invariant relative Lorentz 

factor and 7c is the Lorentz factor associated with velocity of the center of momentum (CM) frame 
relative to the lab frame (Dermer 1984, 1985). Here, u = cos ^u, and the variable ^ is defined as 
the angle between the velocity of the CM frame and the particle 1 velocity in that frame. The 
Kronecker delta 812 prevents double counting. 

To arrive at the relativistically correct reaction rate for two mono-energetic distributions, we 
take the particle functions /*(7i), z = 1, 2, to be 

f:{l^) = ^5H-l^), (11) 

where (ii = (1 — 1/7^)^/^, 7^' are the laboratory frame Lorentz factors of the colliding particles 
expressed in terms of the relative Lorentz factor 7^, and the CM Lorentz factor 7c: 

?ni(2) 7l(2) = ^^17^ [("^2(l)7r + ^-1(2)) ± /3c?7l2(l)/?r7r , (12) 

where S = mf -I- ml + 2mim2 7r and -I- is for i = I, while — sign is for i = 2. 

In deriving the energy exchange coefficient, we use inside the integral of Equation (10) the 
energy exchange rate ((7(7,., 7c)Ai?) for Coulomb interactions averaged over all small (see below) 
scattering angles in the center of momentum frame given explicitly by 

{a{^r,lc)AE) = f dV* ^ AE{f*) , (13) 
J dp* 



- 10 - 



where AE is the energy exchange per scattering (in the laboratory frame) and p* is the particle 
momentum in the center of momentum frame. 

Before we proceed with the derivation, we note that it only makes sense to talk about FP 
coefficients for the Coulomb interactions for small scattering angles. Even though such interactions 
dominate in the Coulomb scattering cross section, one must explicitly exclude the large scattering 
angles in the integrals for the FP coefficients. We will show later that inclusion of these angles 
in the FP coefficients (as done by DL89, for example) leads to erroneous results. Instead, the 
large angle scatterings must be included by means of the exact Boltzmann equation, that is the 
original integral equation from which the FP equation is derived. Following the standard practice, 
we neglect the large scattering angles in the Coulomb interactions in this paper. 

Using the notation of Dermer (1984, 1985), the change in electron energy is 

A£' = 7c/3cP* (cos** - 1) cos /X- sin \E'* cos (;;!>* sin (U , (14) 

where ^* is the scattering angle in the center of momentum system, and (j)* is the polar angle in the 
C-system of coordinates (i.e., the CM frame in which the z-axis is directed along the unscattered 
particle 1 momentum; see Dermer 1985, Fig.l). Because the CM frame Lorentz 7-factor 7* of the 
colliding electrons is 7* = a/ {jr + l)/2 for particles of equal mass, and p* = /?* 7*, the integration 
in Equation (13) is only carried out over angles ^* from to 27r and ** from xE'^jn to some ^E'max 
(see below). 

Since we limit ourselves to the small scattering angles, ^ 1, we can leave only the leading 
terms in the Coulomb scattering cross section, which is then equal to 

da 27rr,2 (27*2 _ 1)2 1 

(15) 



dcOS** ry*2^j*2_l^2 (1-C0S2^'*)2 

(see, e.g., Jauch & Rohrlich 1980), where is the classical electron radius. 

The exact value of usually considered to be unimportant and therefore often chosen as 

7r/2 for simplicity (e.g., Landau and Lifshitz 1981; Dermer 1985). Because the second term inside 
the brackets in Equation (14) gives zero after averaging over cp* , the energy exchange rate in the 
CM frame is then 

{a{jr,lc)^E) = dcos**- — 7c/3cP* (cos**- l)cosM . (16) 

Jcos**. dcosW* L J 

min 

Taking the integral and expressing p* in terms of 7c and 7^-, we obtain 

{aijr,lc)AE) = STrre^ InA ^f^f = u , (17) 



where InA = In 



(l-cos*^^)/(l-cos*;;J 



1/2 



is the Coulomb logarithm. We set InA = 20 

throughout the paper. 

From Equation (10), we now infer the general expression for the energy exchange coefficient 
(note that this expression is now for a single electron, so n2 no longer enters the expression; also 



the factor with the Kronecker ^12 was discarded since we now are interested in the energy exchange 
rate for a given particle rather than the total interaction rate) : 



a(7,7i) 



d7.(7'-l) / d7c(7c-l)'/' / ciu,5{7-72}H7i-7;}(^^(7.,7c)Ai?) . 

(18) 



Performing the third integration first, one finds 



J du6{-f--f!,}{a{jr,1c)AE) 



Sirr^lnA 57^(7, - l)/2 
mi (72- 1)2(7. -l)7c/3c 



7 -7c 



7r + 1 
(-5V2/me) 



(19) 



subject to the condition that | cos/x| = \u\ < 1, which leads to the constraint 



(g^/yme) 



7 -7c 



7r + 1 



< 1 . 



(20) 



According to Equation (12), j^. is fixed by the expression 7c = (7 + li) / {S^^"^ /m^) when mi = 
m2 = rng. Equation (20) then leads to a restriction on the range of physically accessible values of 
7r in the integral of Equation (18), which we write as 7" < 7r < 7"*", where 7^ = 717(1 ± PiP)- 
The next step is to compute the integral over 7c: 



f f Airr^ In A 

J c^7c7c/3c J duSij - 7^} (5{7i - 7;} ((7(7^, jc)^E) = - -^-^—^^j^^——^ (7 - 71) 

1-1 ^ 
The final integral leads to the energy exchange coefficient, which we write as 

2 7rre'^cni In A 



(21) 



where x(7)7i) is the integral function defined by 



(7-7i)x(7,7i) , 



(22) 



7 

X(7,7i) = J 



dx 



X 



^{x + l){x-lf 



(23) 



This integration can be carried out analytically, giving 
X(7,7i) = 



x + l 



X- 1 



+ 2 sinh 



X — 1 



+ 7x2-1 



(24) 



Note that since 7^ is symmetric in 7 and 71, the energy exchange rate given by Eq.(22) is 

antisymmetric with respect to the interchange 7 ^ 71 . 

For the diff'usion coefficient, we proceed in a similar fashion. Now, however, it is imperative 
(see Appendix A. 4) that we choose ^^ax ^ 1- We square Equation (14) and consider the term with 
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(1— cos^*)^. If one takes ^I^max = angles ~ 1 will dominate the contribution of this term to 
the integral over ^'*. But this is the parameter space where the FP equation becomes invalid, thus 
we cannot simply suggest ^'^ax = ''^ 1'^ must use ^^ax ^ 1- ^&\> case, contribution of the 
first term in Equation (14) is negligible compared to the second term in the brackets. Accordingly, 



(c7(7r,7c)(Ai?)^ 



47rreV(7c'-l) 



7.2 _i 

and after carrying out the necessary integrations, we arrive at 



lnA(l 



u 



A-KTe^crii In A 



^i7i'/37' 
The function C(7)7i) is defined as 



^ (7 - 7i)^x(7,7i) + C(7,7i) 



(25) 



(26) 



C(7,7i) 



dx 



(7 + 7i)' 
2(1 + a;) 



(27) 



The integrand is a relatively well behaved function, and the integral is easy to evaluate numerically 
as well as analytically. 

Figures 1 and 2 show the FP energy exchange 0(7,71) and dispersion D(7, 71) coefficients, 
respectively, as functions of 7 for three different values of 71. In these figures, tc is the Thomson 
time divided by the Coulomb logarithm {tc = l/ngCciT In A). It is evident that 0(7,71) is a 
discontinuous function of 7; at 7 = 71 it changes sign from a non-zero positive value to a non-zero 
negative one. This non-physical behavior is caused by the divergence of the original Coulomb 
scattering cross section for \vi — V2\ <^ Vi. This discontinuity in 0(7,71) for |i;i — 'U2I ^ '^i is 
not a problem physically, because the parameter space for such interactions is very small and will 
contribute negligible to the electron-electron Coulomb interactions. It is also never a problem 
numerically. The convolution of 0(7,71) over any real physical distribution (i.e., a reasonably 
smooth one) smoothes out this non-physical behavior and practically eliminates the problem (one 
can show that electrons with energies I71 — 7I = A£' have 0(7,71) oc (AE')^ for any continuous 
/(7))- 

A well known fact is seen from Figure 1: particles with lower energy both gain and lose energy 
(depending on how energetic their interacting partners are) faster than the more energetic ones. 
This happens because of the dependence of the Coulomb scattering cross section on energy of the 
colliding particles in the CM frame (Eq. 14). 

It is possible to derive analytical limits for the FP coefficients. For f3i <^ 1 one gets from 
Equation (22) the limiting form for 0(7,71): 



r (3/2) l/tc/3i, for 7<7i 
«(7,7i) = { 

i -(3/2) l/tcA for 7 > 71 



(28) 



In the limit 7 S> 71 3> 1 it is also possible to have an analytic approximation for 0(7,71): 

0(7,71) ^ -(3/2) (1 - i) ^ -(3/2) 1 1 . (29) 
In V'Yi 'Y/ In Ti 
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We can merge these two expressions in the case of = 7 — 1 3> £"1 = 71 — 1 to give 



(30) 



4.2. FP Coefficients For Arbitrary Particle Distributions: Examples 

Now that we have the mono-energetic coefficients in hand, finding those for any given arbitrary 
particle distribution /(71) is simply a matter of convolving these with the particle profile: 



and similarly for -D(7). 

To illustrate how this works in practice, we show in Figures 3 and 4 the e-e Coulomb energy 
exchange and diffusion coefficients for three different distributions: (1) a perfect Maxwellian at 
temperature ©e = 1 > (2) a Gaussian with the same average energy, and (3) a power-law 7(7) = 
/37~^ (with the maximum 7max — 1 = 100, 7miii — 1 = 0.007 and p = 2.48). The power-law and the 
Gaussian are chosen as two opposite representative cases, the former broader than and the latter 
narrower than the Maxwellian. All three distributions are normalized to unity and are plotted in 
Figure 5. 

We note that the zeros of the energy exchange coefficient occur not at 7 = (71), where (71) is the 
average Lorentz factor of the distribution, but rather at the maximum of the number distribution 
function. It is evident also that the quantity / d'y\a{'y)\ is larger for the more diffuse distributions. 
The diffusion coefficient behaves in the opposite way. The larger diffusion coefficient corresponds 
to the sharper distribution function. Summarizing, we conclude that if the distribution is broader 
than a Maxwellian, it will be brought to the thermal distribution by the energy exchange coefficient, 
whereas if it is narrower, it will reach the thermal distribution as a result of the diffusion coefficient. 

Dermer and Liang (1989) derived the energy exchange and the diffusion coefficients for the 
Coulomb interactions of a test electron with a background Maxwell-Boltzmann electron plasma. 
For the Maxwellian distribution, our results for the energy exchange coefficient are in complete 
agreement (to 1/lnA precision) with those of DL89. For the diffusion coefficient, however, our 
expression disagrees with theirs at both the high and low energy end. Comparing the two 
expressions (Boettcher, 1996, private communication) one finds that they agree with each other well 
in the region E Of, (up to the same precision of 1/ In A) , but seriously disagree for E <^@f, and 
E » Be- We find that r>(7, 0e) const 20e|a(oo, Be)!) for E » 306, while DL89 found that 
Z)(7, 6e) — const' X £^ in that limit. We also find -D(7, 0e) oc on the low energy end, whereas 
DL89 found -D(7, Qe) ^ const. The difference results from the fact that DL89 used ^^ax = 
and also left the terms leading to large scattering angles in the Coulomb cross section (see Appendix 
A. 4) Taking into account the fact that 0(7, @e) — —const for E ^ 36e, we see from Equation (8) 
that the equilibrium Maxwellian distribution can be reached only for D{'y, 6e) — const at the high 
energy end. If instead, -D(7, &e) — const' x E were the correct behavior, then one would find that 
the equilibrium f{j) is a power-law in the high energy limit. Since the original Boltzmann collision 
integral for Coulomb interactions leads to a Maxwellian electron distribution, this property must 
be preserved by the FP equation. 




(31) 
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4.3 The Thermalization Process and Time Scales 

An interesting question to consider is the following: Is self-consistent thermalization faster or 
slower than the thermalization of a test particle distribution interacting with a (fixed) background 
Maxwell-Boltzmann plasma? The answer is not obvious because both of the FP coefficients 
play a role in establishing the Maxwellian distribution, and as we saw in the previous section, 
if a distribution's energy exchange coefficient is larger than that of the Maxwellian, the diffusion 
coefficient is smaller, and vice versa. 

To answer this question, at least illustratively, we consider the thermalization of a Gaussian 
electron distribution in two special cases: (a) when the electrons interact with each other (i.e., 
solving the fully non-linear problem), and (b) when they thermalize on a background MB gas 
with the same average energy as that of the Gaussian distribution. In order to characterize the 
interaction quantitatively, we define the thermalization time scale tr to be the time it takes the 
particle population to reach a deviation of 5% or less from a perfect Maxwellian, where the deviation 
£ is defined to be 

oo 

/ dj{j-l)\f{j,t)-fM{7)\ 

e = 5S • (32) 

/ ^7(7- 1)/m(7) 
1 

In this expression, /{jji) and fuil) are, respectively, the actual time-dependent distribution 
function and the perfect Maxwellian with the same average energy and number of particles. 

Our tests indicate that to within ^ 5% accuracy, the time scale for the test particle 
thermalization coincides with that for the self-interacting case. The explanation for this appears 
to be that most of the elapsed time occurs during the period when the distribution is close to 
Maxwellian, when in fact the FP coefficients for the time-dependent distribution and those for the 
Maxwellian are very similar. As such, our results are in very good agreement with the findings of 
DL89: 

= 4tTep(7ri/2 _ 1.201/4 ^ 20y2)/inA , (33) 

where tx = {nicaT)~^ is the Thomson time. At the same time wc find that the thermalization 
time depends on the shape of the initial distribution. However, if one keeps the shape of the 
initial distribution fixed (e.g., a power-law) and observes the temperature dependence of the 
thermalization time, it is practically identical to Equation (33). 

Figures 6 and 7 show the time evolution of the initial distribution functions, a power-law and a 
Gaussian, respectively, as they approach the equilibrium Maxwellian. The equilibrium temperature 
is the same for both distributions and equals 0.3 mgC^. Figure 8 shows the time dependence of 
the deviation e for both cases. It is notable that thermalization of the Gaussian happens much 
faster than that of the power-law. (When comparing the figures for these two cases, one should 
bear in mind the fact that different parts of the distribution function contribute to the integral in 
Eq. (32) with different weights, namely, E"^ . Thus, the high energy end of the electron distribution 
contributes much more to the deviation e than the low energy end.) To explain this fact, we 
should consider the Coulomb interaction time scale for individual particles. Let us define the 
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K7)I + ^J, (34) 



particle energy time scale as 

where 0(7) and .0(7) are the Maxwellian FP coefficients for the equilibrium electron temperature. 
The physical sense of tg is that it is the time it takes to change the particle energy considerably 
(say, by a factor of 2). To obtain the thermalization time scale for a distribution as a whole, one 
should average l/te(7) over the distribution and take the inverse of this value. Applying this to 
the two thermalization processes considered (Figs. 6 & 7), the power-law has substantially more 
particles at the higher energy end than the Gaussian, and it therefore has a longer thermalization 
time, in accordance with Figure 8. 

Figure 8 also shows that the equilibrium distribution function does not exactly coincide with 
the Maxwellian since e does not go to zero when t — > 00, but rather stays constant at ~ 10~^. 
This occurs because the diflFerencing scheme (see Eq. 2) , modifies the true FP coefficients in order 
to preserve the energy and number of particles (see Appendix). One can introduce more elaborate 
energy-dependent corrections to the FP coefficients, and it will likely reduce the deviation e by an 
order of magnitude or more. However, since a precision of about 0.01 is satisfactory for our goals, 
we will leave these improvements to the future. 

5. Collision Integrals For Electron-Proton and Electron-Photon Interactions 

5. 1 . Fokker-Planck Coefficients For Electron-Proton Interactions 

The fact that electrons and protons are subject to the same Coulomb force as electrons with 
electrons and electrons with positrons (other than for the cross section) allows us to use the same 
Fokker-Planck approach as described above to account for these interactions. The FP coefficients 
for electron-proton (e-p hereafter) scatterings will differ from those of lepton-lepton interactions for 
the additional reason that the particle masses are no longer the same. In this section, we shall derive 
the e-p energy exchange ap{'y,Ep) and diffusion Dp{'j,Ep) coefficients along the lines established 
above (Ep is the proton kinetic energy in units of nipC^). Unlike earlier work in which the FP 
coefficients were derived for a thermal proton distribution (DL89), we here solve the problem for 
an arbitrary proton proffie. 

We follow the same strategy as above and neglect terms of order lower than In A in the final 
expression. Assuming infinitely massive scatterers (the protons), the e-p differential cross section 
in the center of momentum frame is 

dcos^'* 2 p^-f^sin'^{^*/2) ^ ' 

(see Jauch & Rohrlich 1980; we neglect the factor (/3r sin ^I'*)^ compared to unity in the numerator, 
since we again only include ^ 1). The validity of this treatment is limited to 7^, ^ nip/ 
where nip is the proton mass. The mono-energetic e-p energy exchange coefficient is found to be 

«p(7,7p) = p^2fj ^2 ] (^2 _ 1)3/2 ~ ^){mp7p - me7) + [mp + nie){jp - 7)] , (36) 
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where 7^ and 7 are the proton and electron Lorentz factors, respectively, 7"*= = 77p(l ± P(3p), and 
Sp = nip + ml + 2mpme^r- We note that due to the difference in electron and proton masses, 
ap(7, 7p) 7^ even when 7 = 7^. The energy exchange coefficient ap(7, 7p) is shown in Figure 9 
as a function of electron energy, for three different values of the proton energy. Note that there is 
still a discontinuity at 7 = 7^, however. 

The corresponding e-p diffusion coefficient is given by 



2Trrlcnp In Anipnie f d^rlr^p 



(mp7p + me7)^ _ [(7r - l)(fflp7p - ffl-e7) + ("^p + "^e)(7p - 7)]^ 

Sp Spin^^-i) 



(37) 



To get the c-p FP coefficients for an arbitrary proton distribution fpijp), one integrates 
Equations (36) and (37) over fpi 

00 

apil) = J d'yp ap(7, 7^) fpijp) . (38) 
1 

We plot the e-p diffusion coefficient for Maxwellian protons with a temperature kTp = 10 MeV and 

100 MeV in Figure 10. For comparison, wc also show the electron-electron diffusion coefficient for 
a Maxwellian distribution with kTf, = 1 rrieC^. A graph of the the e-p energy exchange coefficient 
can be found in DL89. Note that the e-p FP coefficients differ from those of the electron-electron 
interaction roughly by a factor {me/mp)Tp/Tf,. 

5.2. Electron-Photon Interactions 

In this section, we present the electron-photon collision integral corresponding to Compton 
scattering. This part of the kinetic equation has already been dealt with extensively in the literature 
(see references below in this section), especially for isotropic photon and particle distributions. 
With the assumption of isotropy in both the electron and photon spectra, the electron-photon part 
of the collision integral is 

i^^^)e, = "^^^^ / ^^^M^(^'7) + // daj'dj'R{u',j')P{r,7',u;')N{u')f{j') , (39) 

where R{uj,^) is the scattering rate between photons of energy uj and electrons of energy 7; 
P(7;7',cj') is the probability of scattering an electron with initial energy 7' to a final energy 
7 in a collision with a photon of energy u)', and N[u)) duj is the number of photons per unit volume 
with energy between uj and uj + duj. For the angle-averaged scattering rate R{aj,^), we use the 
expression given in Equation (2.3) in Coppi & Blandford (1990). 

The probability P(7;7',u;') was first derived by Jones (1968) and later corrected by Coppi & 
Blandford (1990). We carried out the last integral in the expression given by Coppi & Blandford 
(1990) to get P(7;7',a;') integrated over all angles. We do not give the rather lengthy expression 
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since more complete work exists in the literature (e.g., Nagirner and Poutanen 1994 and references 
therein) . We checked the validity of our expression by comparing it with the numerical integrations 
of Coppi & Blandford (1990) and by integrating P(7; 7', oj') over all the accessible values of 7, which 
should give unity. 

A problem that may arise when Equation (39) is integrated numerically is that the electron 
energy bin size is finite, which can lead to a situation where scatterings in the low energy regime 
cannot correctly transfer the particles in energy space. Consider a photon with energy oj (in units 
of rrieC^) scattering in the Thomson limit, so that 0)7^ <^ 1. The photon emerges with energy 
uj' ~ 7^0; after one scattering oflF an electron with a Lorentz factor 7. If the change Acj = w' — uj'm 
the photon energy is smaller than the electron energy bin size, such a scattering will 'fall through' 
the energy space mesh. When the density of low energy photons is large, the error introduced by 
this effect can be large. 

Because of limited computer resources, it is often impractical to simply reduce the energy bin 
size. Instead, the photon energy range used in Equation (39) can be divided into 2 sections. Let 
us define a break energy a;;, (7) in the photon spectrum by the expression 0Ji,{'y) = min{^ (7 — 
1), 3/(47)}. For low energy photons, that is for uj < uj^, we introduce the alternative cooling rate 

jc{l)= J dioN{iu)Riiu,-f){Mj))-io) , (40) 


where {tOsij)) is the mean scattered energy of the photon lo interacting with an electron with 
Lorentz factor 7. This cooling rate is in fact the electron-photon energy exchange rate and is 
analogous to the quantities 0(7,71) and ap(7,7p). 

One must also take into account the electron diffusion associated with Compton scatterings, 
that is the full diffusion coefficient must include a contribution from electron-photon scatterings 
as well. We define the electron-photon diffusion coefficient as 

Dc{l)= J dujN{cv)R{cv,j) {{cvl{j))-{u;,{j)f) , (41) 


where {ijJsi^)) is the second moment of the distribution of scattered photons in the scattering of 
the photon u and electron 7. 

We use the full kinetic equation when the photon energy is a; > ^(,(7). In this case, there is 
no problem with the finite energy bin size because the change in photon energy when these higher 
frequency photons scatter is large and A7 (7 — 1). Equation (39) is then written in the form 



\ dt J 



67 



2 ^max 

--^[icf{l,t)] + ^^[Dci7)f{7,t)]-f{j) I du:N{u)R{u,^) 



^max 

+ y duj' Jdj'Riu;',j')P{r,7',u;')N{u;')f{j'), (42) 



1^6 (7) 



where comax is the maximum photon energy in the spectrum. To compute both 7c (7) and Dc (7) 
numerically we use a much finer resolution in photon energy space than that for u > uJbi'j)- 

We have tested the validity of this approach by first varying the 'break frequency' uJb{'y) and 
found that the variations in the results are always much smaller than the chosen relative variation 
in cuh{'y). We also varied the number of bins N for a fixed electron energy range and found that 
this too makes a negligible difference on the results. These tests would unambiguously point to 
inconsistencies if they were there. 

Finally, for pair annihilation we use the rate i?(7_,7+) given by Svensson (1982), where 7_ 
and 7_|_ are electron and positron gamma-factors, respectively. For pair creation we adopt the 
expressions in Coppi &; Blandford (1990). 

6. Tests 
6.1. The Fall Equation 

Our goal in this section is to demonstrate the use of the self-consistent FP approach in 
the treatment of photon-particle interactions in astrophysical plasmas. Situations where one 
can assume that a part of the distribution is Maxwellian have been studied in some detail for 
the case of non-thermal models (e.g, Zdziarski et al. 1990; Svensson 1994, and references cited 
therein). GhiscUini, Guilbert & Svensson (1988) studied the thcrmalization of the low energy 
end of the non-thermal distribution by synchrotron reabsorption; they neglected e-e Coulomb and 
inverse Compton interactions. DL89 developed a perturbative approach to the FP equation for 
e-e Coulomb interactions, but have not applied the method to problems that require finding the 
photon distribution self-consistently. 

We first present the complete equation, including all the processes discussed above, i.e., 
Coulomb and Compton interactions, as well as pair creation and annihilation. Wc shall not 
here attempt to make a detailed study of the various applications, which would require a careful 
examination of source geometries and a thorough search in parameter space. This subject is vast 
and will be covered in future publications. 

All the distributions are assumed to be spatially uniform inside of a spherical volume with 
radius R. Since the number of positrons is always different (i.e., smaller) than the number of 
electrons, the positron distribution is never identical to that of the electrons. This is true even 
though the electron-electron and electron-positron cross sections are the same within the scope 
of the approximation we use here (see §4.1). The processes that shape these particle profiles 
differently include pair annihilation and, when included, bremsstrahlung emission and pair escape. 
In the following, we neglect the contribution of bremsstrahlung to the overall photon production, 
assuming instead that the soft photons are injected by an external source. Index 1 will be used to 
denote electrons, whereas 2 pertains to positrons. 

The full equation for the electrons reads 



+ 



67 



i?i(7)/i(7,t)+ 5(7) 



(43) 



- 19 - 



where ^(7) is the total energy exchange rate for an electron with Lorentz factor 7; that is, ^1(7) 
is the sum of energy exchange rates from processes which can be treated using a Fokker-Planck 
approach. This excludes inverse Compton interactions which are instead incorporated into the 
term (5/1(7, i)/5t)e7, given in §4.2. Similarly, .0(7) is the electron diffusion coefficient, a sum 
of the electron-electron, positron-electron and proton-electron diffusion coefficients. 
The function Ri (7) is the electron annihilation rate, 

^1(7)= J d-i2f2{l2)R{l,l2) , (44) 

1 

where i?(7,72) is the angle-averaged annihilation rate for an electron with energy 7 and positron 
with energy 72 as given by Svensson (1982). 5(7) is a 'source' term that represents pair production 
by photon-photon collisions. The pair creation rate is an integral over the photon distribution 
N{u,t) as given by Coppi & Blandford (1990). 

The full equation for the positrons has an identical form, except that the subscripts 1 and 2 
are interchanged. The often used quantitative measure of the heating power going into electrons 
is the dimensionless 'hard' compactness parameter 

lh = ^^^^ j d^A+{^,t) [/i(7,t) + /2(7,i)], (45) 

where ^'''(7, t) is the energy and time dependent heating rate for an electron/positron with energy 
7. For the proton heating rate given by Equation (38), the time dependence results from the 
change of the proton temperature in time. Letting Ih represent the energy transfer rate to the 
protons (due to dissipational processes), the proton temperature is then found by balancing this 
heating and the Coulomb cooling by cold electrons at every instant in time. 
The equation for the photons reads 

^ = iV(.) - (^) + P(.) - iV{a,)/t.„(a,) , (46) 

where N{u)) is the rate of external photon injection (photons are here assumed to be injected 
uniformly in space throughout the plasma containing region), P(w) is the number of photons with 
energy oj created per unit time per unit volume per unit energy by pair annihilation, and tesc{'^) 
is the photon escape time as given by Lightman &; Zdziarski 1987: 

= - [1 + (1/3) (t>{u:) TKN(a;)] , (47) 
c 

in terms of the size R of the region. The energy dependent optical thickness rKN(w) is here defined 

by 

t{u) = TTCrKN(t^)/crT , 

where ty,^{lo) is the Klein-Nishina cross-section (see, e.g., Jauch and Rohrlich 1980, Eq.[11.24]). 
The function 0(aj) takes into account the fact that forward scattering for high energy photons 
becomes important (see Lightman and Zdziarski 1987; Equation [21b]). We should mention that 
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this treatment of the photon radiative transport is only approximate, especially for a low optical 
depth. However, the solution found this way will undoubtedly bear the same qualitative imprints 
of the 'exact' electron distribution as the still to be found solution which includes both radiative 
transport for the photons and the FP treatment of the electrons. 

The term [dN{LL!)/dt]e'y represents the full Compton interaction with electrons and positrons 
and is given by the collision integral 

OO OO 

-(^^)^^ = -^M JduJ d7i?(cu,7)[/i(7) + /2(7)] 

1 

OO OO 

+ Jdu;' J dj'P{j,j',u')Riu;',j')Niu;') [Mj') + h^')] , (48) 

1 

where R{uj, 7) and P(7, 7', t^^') are defined in §4.2, and 7 + = 7' + to' . 

In the tests reported here, we take the photon injection rate A'^(ti;) to be that of a blackbody 
spectrum with dimensionless temperature Qb = kTb/rUeC^ in the range = 10~^ — 10~*. This 
injection rate is characterized by the dimensionless 'soft' photon compactness parameter (e.g., 
Lightman k, Zdziarski 1987) 



'~ 3c 



■jduoujN{u). (49) 



We also need to specify the electron heating mechanism. It is well known (Guilbert, Fabian k, 
Stepney 1982), that the electron-proton Coulomb interactions are very inefficient, and if cooling is 
effective, it is impossible to account for the observed high electron temperature 9e ^ 0.2 observed 
in Seyfert Galaxies, as implied by their high exponential cutoff energy (e.g., Svensson 1996 and 
references therein). One then needs to introduce some ad hoc heating mechanism for the electrons. 
Many workers do not specify the mechanism by which the energy is transferred to the electrons 
and instead specify the heating rate per electron (independent of the electron energy) per unit 
volume. This is indeed the only information needed if one assumes that the particles are thermal. 

However, for our purposes we need to know the diflfusion coefficient for this heating interaction. 
Since it is beyond the scope of this paper to explore the full range of interactions that heat 
the electrons in real astrophysical plasmas, we model the electron heating by the e-p Coulomb 
interactions with a fixed proton temperature, kT^ = 20 MeV, but in addition, we multiply this 
interaction (both the energy exchange and diffusion) by a factor (equal to Tp/To, see § 6.4.c) which 
is needed to account for the electron heating rate given by the specified Ih- While this is obviously 
a non-physical interaction, it will allow us to make some progress in understanding what one can 
expect to happen with the electron distribution. It also will allow us to have a 'reference scale'; 
whether the electron distribution is actually broader or narrower than the one found in this test 
model depends on whether the actual electron heating mechanism produces comparatively more 
or less diflFusion. 

6.3. Results 

In this section we apply our technique to a test situation where the electrons and positrons 
are heated by some external source as described in section 6.1, and cooled by inverse Compton 
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scatterings with soft intense external radiation. The description of the numerical procedure used 
to solve for the exact electron distribution can be found in Appendix A. 3. Ghisellini, Haardt 
and Fabian (1993) (GHF hereafter) have shown that the exact shape of the electron distribution 
produces very little effect on the photon spectra as long as the multiple Compton scatterings are 
the main emission mechanism. In particular, they considered three different particle distributions 
that had the same normalization and the same value of (7^ — 1): a Maxwellian, a power-law, and 
a constant (i.e., 7(7) =const). The resulting spectra were hardly distinguishable. 

We conduct similar tests here. For the test electron distribution, we choose: 1) the exact 
electron distribution found by using the techniques developed in this paper; 2) a Maxwellian; 3) a 
distribution broader than a Maxwellian (see Figs. 11a and 12a). In the last case, the functional 
shape of the distribution function is given by f{E) = const a/ fmxw{E, &e), where finxw{E,Qe) is 
a Maxwellian distribution with temperature Q^. Just like in the familiar Maxwellian case, there 
are again two parameters, ©e and the normalization, which are found self-consistently. One very 
important difference of our tests from those of GHF is that we find the equilibrium pair number 
density by solving the exact pair balance, so that the relative positron number density, z = U-^-fup, 
is different for the three tests. This results in a different optical depth for the different distributions, 
even though the proton optical depth Tp is the same for all three. 

We first discuss the photon spectra. Figures lib and 12b show the photon spectra for Ih = 
Is = 420, Tp = 0.05 and Ih = 8.4, Ig = 2.1, Tp = 0.02 respectively. The input and output parameters 
of these tests are given in Table 1. For the case of the higher optical depth (Fig. lib), it is readily 
seen that the spectrum does not vary much from one case to another, at least in the X-ray energy 
range (2-18 keV). This is similar to the findings of GHF, and is explained by the fact that the 
spectrum is formed by many repeated scatterings, which simply form a power-law after two to 
three interactions. 

However, for the second test our results disagree with conclusions of GHF. We find qualitative 
differences between the Comptonized spectra produced by the three test electron distributions 
(Fig. 12b). The difference with GHF is caused by the fact that the calculated here exact electron 
distribution is narrower than a thermal one (Fig. 12a), and also that we test the parameter space 
with a lower optical depth. GHF considered only a subset of the parameter space available for 
the distribution function; they only considered cases where the electron distribution is broader, 
more diffuse than the thermal distribution. At the present, it is completely unclear if this is indeed 
what happens. The answer to the question regarding the electron distribution function can only 
be given if one specifies the interactions affecting the electrons. 

6.4. The Strongly Non-Majcwellian Electron Distribution 
a. The electron distribution profile 

Using the tests described below as examples, we will discuss the parameter space where A is 
much larger than unity, i.e., where strongly non-thermal distributions are expected. Figures 11a 
and 12a show the three test electron distributions for the tests described in 6.3 (A c± 80 and 10 

for the two figures, respectively). The exact electron distributions are shown as solid curves, and 
are much sharper than the thermal electron distribution (dashed line). For the parameter space 
considered here, the electrons with energy E <C (£"), where {E) is the electron average energy, 
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interact mostly with protons which strive to push these electrons to higher energy. The electrons 
with energy E S> {E), similarly interact mostly with photons. Therefore, the electron distribution 
appears to be squeezed by the photon cooling on the high energy end and the proton heating on 
the low energy end (see Fig. 13). 

As surprising as this might sound, the exact distributions in Figures 11a) and 12a) are actually 
as simple to find as the familiar Maxwellian distribution (corresponding to the given compactnesses 
and Tp). The FP part of the full electron Equation (43) is by far the dominant one. That means that 
the actual shape of the particle distributions is controlled by the FP part, while the normalizations 
are governed by the much slower annihilation and creation processes. Thus we can find the shape 
of the distributions solving just the stationary FP equation. Furthermore, for A 3> 1, both electron 
heating and cooling greatly exceed the c-c Coulomb interactions (for an example see Fig. 13). 
Essentially, the FP equation becomes linear in the electron distribution function, which greatly 
simplifies the solution of this equation. The stationary (i.e., with zero time derivative) FP Equation 
(8) has then an exact solution: 



„, , const 
^^^^ = D{^) 



dj' [A{j')/D{j')] 



(50) 



We have dropped subscripts 1 and 2 on the distribution function since in this approximation the 
positron distribution has the same shape as that of the electrons, and only the normalization is 
different.^ 

Figure 14 shows the exact distribution (same as in Fig. 12a) and various fits to it. The heavy 
dashed curve shows the function given by Equation (50) when ^(7) and .0(7) include all the 
interactions. The larger dots show the same Equation (50) but with the e-e interactions excluded, 
i.e., the linear case. As one can see, both fits are very good. And the faint dotted curve gives an 
analytic fit to Equation (50). 

We now explain this fit. Figure 13 shows all the FP coefficients for the test previously shown 
in Fig. 12. Notice that the proton diffusion coefficient is almost a constant, so we can approximate 
D(7) c± .0(70) = Do (where 70 is the point where electron cooling is equal to electron heating). 
We also approximate the electron-proton energy exchange coefficient as a constant, A"*" (7) ~ ^0 = 
y4+(7o), which should be sufficiently accurate around the peak energy Eq (see the Figure). Using 
the fact that ^"(7) oc /?^7^, we can write ^(7) ~ Aq[1 — (/?7//3o7o)^]- Using this expression in 
(50), we get 

2An 



7(7) ~ const exp 



V(7o^-7V3) 



(51) 



Notice that this expression is also quite accurate around the peak of the exact distribution function, 
but deviates on the low energy end. This deviation will be negligible for the distribution as a whole 
if the electron heating is more or less energy independent. 



t This similarity in the electron and positron distributions exists in the exact simulations as 
well, and we therefore do not show the positron curves in Figures 11a &: 12a. 
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One can easily show that the dispersion of the electron distribution, defined as (A7) = [((7 — 
70)^)]^^^, where the averaging is taken over the electron distribution, is 

1/2 

(52) 

Therefore, for the same heating rate (which will lead to the same 70), it is the diffusion coefficient 
that determines the dispersion of the distribution. Physical conditions under which the 'wide' 
electron distribution might form can now be understood in terms of a very large diffusion coefficient 
which might result from some particular form of the electron interaction. 

6.4.b. Normalization of the electron distribution 

The electron-positron annihilation rate depends rather weakly on the energy of the colliding 
particles as long as their energy is smaller than a few MeV (e.g., Coppi and Blandford 1990, 
Equation 3.7). Therefore, the exact particle profile does not change the pair annihilation rate 
substantially. Nevertheless, its importance shows up in the shape of the photon spectra around 
the pair producing energy (see Figs, lib and 12b). We can indirectly find the equilibrium pair 
number density by considering the Compton y-parameter (see, e.g., Rybicki and Lightman 1979). 
It is defined as 

y = {G) (tesc)AT, (53) 

where {G) is the average fractional gain in the photon energy in one scattering, (tesc) is the 
average escape time and tr is the Thomson time. By definition, the electron cooling rate is 
{dE/dt)~ = caT{G)Uph, where Uph is the photon energy density in the source. In equilibrium, the 
total electron heating, L^, equals the total cooling, so that 



(A7) ~ 



Lh = Voripil + 2z) {dE/dt)- = carUpil + 2z){G) Vq Upi,. (54) 

where Vq is the source total volume and z = n^/up, ratio of positron number density to that of 
protons. Using Equation (53) and the fact that L = C/phVo/(tesc)) where L = + is the total 
luminosity of the source, we get 

This derivation did not make any use of the particular shape of the electron distribution, and 

therefore is valid for an arbitrary one, as long as the inverse Compton scattering is the main 
mechanism for the transfer of energy from the particles to the escaping photons. We note that the 
numerical value predicted by Equation (55) should be corrected by taking into account relativistic 
Klein-Nishina decline in the Compton cross section. The principal effect then is to reduce in 
the denominator, since these are the photons that have largest energies. This correction is rather 
trivial to do, and we leave it to the reader. 

Now, since we know the electron distribution shape (Eq. 50), we can find (G) by simple 
convolution of this shape with the particle-photon energy exchange rate. Equations (55) and 
(53) then determine the Thomson optical depth and so the desired normalization of the electron 
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distribution. This is particularly simple for the Thomson scattering limit, e.g., for a spherical 
geometry, 

y ~ tt{1 + rT/3)(4/3)(/3V)- (56) 

For the tests shown in Figure 11, the y-parameter must equal 1/2. Numerically finding (/3^7^) in 
Equation (56), we get y = 0.5 ± 0.01 for all three tests. 

We have shown that, once the electron energy exchange and diffusion coefficients are known, 
it is a simple matter to find the electron distribution. Although we have considered only systems 
where the electron cooling is dominated by the Comptonization of the soft external photons, one 
can extend the methods described here to more complicated situations. We intend to cover this in 
future publications. 

6.4.C. Proton temperature 

The last column in Table 1 shows the equilibrium proton 'temperature'. We define the proton 
temperature such that the e-p heating is equal to {Tp/To) ap(7,To), and fcTo = 20 MeV. For both 
tests, we see that Tp varies from distribution to distribution, and the pattern is such that the more 
diffuse the distribution is, the colder the protons become. 

The explanation for this is rather straightforward. The e-p energy exchange rate (see Fig. 9 in 
this paper and Fig. 10 in DL89) is a decreasing function of the electron energy. It is mostly the low 
energy electrons that are heated by the e-p interactions. Since more diffuse distributions have more 
electrons at low energy, they absorb more energy per electron from the hot protons. Additionally, 
we also saw in the tests we conducted that the more diffuse distributions are associated with the 
larger optical depths. That is, there are more leptons per proton, which again lowers the proton 
temperature. 

Yet another consequence of the change in the electron distribution function is an associated 
change in the proton cooling time scale. Let us define this time scale by the equality Ep/tp = Ap, 
where Ap is the proton cooling rate, and Ep is the proton average energy. Assume that the 
proton heating rate (e.g., from the dissipation of gravitational energy) is a constant. Then Ap is 
independent of the electron distribution, since proton cooling and heating are equal to each other 
in equihbrium. Also, Ep oc Tp, the proton temperature, and thus the proton cooling time scale is 
directly proportional to this temperature. 

Even though we here assume a non-physical form for the proton heating, it is quite likely that 
a similar result will hold for real interactions that transfer energy from the hot protons to the 
electrons in astrophysical plasmas. Indeed, the e-p energy exchange rate must be a decreasing 
function of the electron energy, since after all it must become negative after the electron energy 
surpasses that of the protons. The second part of the argument, based on the change in lepton 
number, holds true irrespective of the exact nature of the e-p interaction. Finally, the proton 
cooling time scale will remain unchanged only when the proton-electron energy transfer rate is 
independent of the proton energy, which is rather unlikely. 

8. Conclusions 

We have presented the general FP method for finding the electron distribution in non- 
Maxwellian hot, isotropic plasmas. The technique can be applied to both time dependent and 
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stationary problems. Approximate expressions (Equations 50 and 51) for the exact electron 
distribution function are developed. These approximate expressions enable one to find the exact 
electron distribution in about the same number of steps as is required to find the temperature and 
equilibrium number density for a Maxwellian electron distribution. 

A simple parameter (§ 2) specifies the section of the parameter space where the treatment 
presented here is important. In general, this occurs for optically thin, mildly relativistic and 
relativistic plasmas that are not too strongly magnetized. Two representative tests show that 
the parameter space is further divided into two regions (the exact boundary between these two 
regions depends on the particle distribution and the as-yet-unknown physics of electron heating; 
we intend to investigate this boundary in future work). In the first region the exact shape of the 
electron distribution does not influence the produced photon spectra, which are close to power- 
laws. While this finding agrees with the work of GHF, we find that the difference in the electron 
profile can lead to different equilibrium pair number densities. This fact will be important in the 
radiation pressure dominated plasmas, since the radiative pull on the protons is proportional to 
the number of leptons per proton (= 1 -|- 2z). Of equal importance is the change in the proton 
temperature. A smaller proton temperature increases the importance of the radiation pressure, 
and for the tests we have presented here, the combined effect is such that the ratio of radiation 
pressure to the proton pressure changes by more than a factor of 10 from the exact to the 'wide' 
distribution. With such a large difference in this ratio, there can be situations where the plasma 
can be either gas or radiation dominated depending on what distribution electrons have. Therefore, 
we conclude that it is probable that the equilibrium (including hydrostatic) plasma state will still 
be sufficiently different for different particle distributions to produce noticeable differences in the 
photon spectrum; this question needs to be investigated further. 

In the other, generally optically thinner, region, the photon spectra bear strong imprints of 
the exact shape of the particle distribution, and are not simple power-laws. The local X-ray index 
and break energies depend on the exact particle profile. Use of the exact electron distribution can 
then limit the parameter space for a particular astrophysical model differently than the use of the 
thermal distribution does. 

A shortcoming of this paper is the use of the diffusive escape time approximation (Equations 
46, 47) instead of the exact radiation transport approach (e.g., Poutanen and Svensson 1996). One 
can, however, include radiation transport in the same way as is done for Maxwellian electrons. 
We anticipate that the 'non-power-law' features in the photon spectra may be even stronger with 
the inclusion of radiation transport. Our treatment of radiation employed the angle averaged 
corrected Jones expression for the Compton scattering. If one is fixing the viewing angle instead, 
the amount of dispersion in the photon upscattered energies become smaller (because one angle 
average is removed) , and the spectra will track the electron distribution better. 

Another interesting question, which should be addressed in future work, is whether there are 
generic limits on the ratio of the electron heating diffusion coefficient to the energy exchange one, 
independent of the exact interaction which heats the electrons. If this is so, then we should be 
able to determine just how much the exact electron distribution can be broader or narrower than 
a Maxwellian. 
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Figure Captions 

Figure 1. Coulomb energy exchange coefScient 0(7,71) for a test particle scattering off a mono- 
energetic electron distribution (i.e., /(7) = d[y — 71]) for three diflferent values of 71. The curves 
are labeled by their corresponding values of 71. The Coulomb logarithm here and throughout the 
paper is set to 20. The time tc is defined as tc = {necax In A)~^, and rie is the total lepton number 
density. 

Figure 2. Coulomb diffusion coefficient D(7,7i) for a test particle scattering off a mono-energetic 
electron distribution (i.e., 7(7) = S^y — 71]) for three different values of 71. The curves are labeled 
by their corresponding values of 71 . 

Figure 3. The electron energy exchange coefficients 0(7) for three different electron distributions 
(plotted in Fig. 5 below). The solid curve corresponds to a Maxwellian distribution with 
temperature &e = kTe/nieC^ = 1. The distributions have the same average energy and are 
normalized to unity. The dotted and dashed curves correspond to a power-law (/(7) = l3j~^; 
p = 2.48) and a Gaussian, respectively. 

Figure 4. The electron diflfusion coefScient .0(7) for the same distributions as in Figure 3. 

Figure 5. The three electron distributions used to compute the FP coefficients plotted in Figures 
3 & 4. All three functions are normalized to unity and have the same average energy. Solid, dotted 
and dashed lines correspond to Maxwellian, power-law and Gaussian distributions, respectively. 
The temperature of the Maxwellian is ©e = 1- 

Figure 6. Time evolution of an initial power-law (/(7) = 7"^; p = 3.72) distribution of electrons 
coming into equilibrium (i.e., approaching a Maxwellian distribution) via Coulomb interactions 
with itself. The temperature of the equilibrium Maxwellian is 0.3. The box on the left shows the 
time corresponding to each curve {tc is defined in the caption of Fig. 1). The solid curve with 
t = 00 corresponds to the perfect Maxwellian. Note that our choice of times to plot f{'-f,t) was 
dictated by convenience of viewing and not by equal spacing in time. In fact, thermalization is 
always slower at later times. 

Figure 7. The same as Figure 6 except that the initial distribution is a Gaussian with the mean 
energy corresponding to a Maxwellian temperature of 0.3 rrieC^. 

Figure 8. Time evolution of £ (defined in Eq. 32), the deviation from a perfect Maxwellian. 
Solid and dotted lines correspond to thermalization of the power-law (Fig. 6) and Gaussian (Fig. 
7). Note the Gaussian relaxes to the thermal distribution much faster than the power-law. For 
comparison, the electron-electron thermalization time scale given by Equation (33) is 0.06 tx for 
the given equilibrium temperature of 0.3. 
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Figure 9. The mono-energetic electron-proton energy exchange coefficient ap(7,7p) as a function 
of the electron energy. The curves are labeled by their corresponding values of the proton kinetic 
energy. 

Figure 10. The electron-proton diffusion coefficient for the case of Maxwellian protons. The 
Maxwellian electron-electron diffusion coefficient for 0e = 1 is shown with a dashed line for 
comparison (same as the solid curve in Fig. 4). 

Figure 11. a) The equilibrium electron distributions for the tests described in section 6.3. The 
parameters are: 1^ = Ig = 420, Tp = 0.05 (also see Table 1). The exact distribution is the 
solution of the full kinetic equation, which includes electron heating, cooling and pair creation and 
annihilation. The Maxwellian curve is the solution of the same problem but assuming that the 
distribution is thermal. The 'wide' electron distribution is again a solution of the same problem, 
but assuming a wider than thermal profile for the distribution. The three distributions have a 
different 'temperature' and number of particles, but produce almost identical spectra (Fig. b). b) 
The photon spectra corresponding to the electron distributions plotted in Fig. 11a. 

Figure 12. a) The same as Fig. 11a, but now for Ih = 8.4, 1^ = 2.1 , Tp = 0.02. b) The 
corresponding photon spectra. Notice that the spectra are distinctly different from each other, 
and are not simple power-laws anymore. See Table 1 for additional information about input and 
output parameters. 

Figure 13. a) The electron energy exchange coefficients plotted for the tests shown in Fig. 12. 
The solid, dotted and dashed curves show the electron-photon cooling (the dominant FP part, see § 
4.2), e-p heating and e-e energy exchange coefficients, b) Same as a), but for the electron diffusion 
coefficients. Note that both energy exchange and diffusion are dominated by electron-photon or 
electron-proton interactions. 

Figure 14. Various fits to the exact electron distribution (shown with the solid curve) for the 
test presented in Fig. 12. The dashed curve corresponds to the function given by Equation (50) 
with all the FP coefficients included, and the big dotted curve corresponds to the solution with 
the electron-electron FP coefficients excluded. The 'faint' dashed curve shows the fit given by the 
simple expression in Equation (51). See text for details. 
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Table 1. The output and input parameters for the tests presented in section 6.3 and plotted in 
Figures 11 & 12. 



ModeP 


h 


h 






Trp 




kT^ (MeV) 


exact 


420 


420 


0.05 


10-4 


0.22 


0.58 


2400 


wide 


420 


420 


0.05 


10-4 


0.41 


0.29 


190 


thermal 


420 


420 


0.05 


10-4 


0.27 


0.48 


1100 


exact 


8.4 


2.1 


0.02 


3 X 10-5 


0.14 


1.46 


440 


wide 


8.4 


2.1 


0.02 


3 X 10-5 


0.20 


0.87 


53 


thermal 


8.4 


2.1 


0.02 


3 X 10-5 


0.17 


1.27 


240 



^ The assumed or exactly determined shape of the electron distribution (see § 6.3). The next four 
columns give the input parameters for the various models: hard and soft compactness (Eq. 50 
& 54, respectively); Tp, the 'proton optical depth'; the injected blackbody photon temperature, 

^ Thomson optical depth of the plasma, tt = Tp(l + 2z), where z is the ratio of positron number 
density to that of the protons. 
^ Average electron energy. 

Proton 'temperature', defined in § 6.4. 
Note that differences in the equilibrium optical depth, and the electron and proton 'temperatures' 
occur even when the radiation spectra are almost indistinguishable. 
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Appendix 



A.l The Finite Differencing scheme 

Finding a numerical scheme to solve Equation (8) that is stable and preserves the essential 
physics (i.e., one that conserves the particle number and energy, and that leads to a perfect 
Maxwellian in equilibrium) can be a non-trivial problem. We shall mention that the need to 
conserve the number and energy is not simply an aesthetic goal, but comes rather from a need 
to address practical problems involving pairs, photons and protons (see § 6). When the escape 
rate of particles becomes small (i.e., when the optical depth is large), we find severe problems 
with all differencing schemes that do not conserve the energy and number "exactly" . In particular, 
unphysical pair runaways may ensue. We shall therefore present the energy and particle conserving 
scheme in detail. The scheme we have settled on is 



J) 



n+l 
k 



-if I 



k+Uk + l 



-Jl 



k-1 



Axk-i + Axk 
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where 
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* aAxk{Axk-i + Axk) ' ^^"^^ 

Here a{E) and D(E) are the energy exchange and diffusion FP coefficients, respectively, which are 
given as integrals over the distribution function / itself. 

The index k refers to the energy (xk is a point in energy space), and n denotes the time. 
The parameters a and p have been introduced to make sure the scheme preserves the particle 
normalization and energy. In our calculations, we have chosen a logarithmic energy spacing, and 
one can show in this case that in order to preserve these aspects, one must choose 



a 



1 + 



_2q_ 

1 + 



(A3) 



where q = Xk+i/xk- To demonstrate that this must be the case, we need only consider a numerical 
realization of the integrals / dE [/"(£;) - f'+^E)] and / dEElf'iE) - f"+\E)] for the evolving 
distribution, and require these to give zero. 

Because Equation (8) is a flux equation, a sensible constraint to impose at the boundary is 
a mirror condition at both the low and high energy ends. This means modifying the differencing 
scheme (Eq. Al) for a few boundary points. We explicitly set the value of the distribution 
function in the first and last points to zero. That means that the two extreme points in the actual 
distribution function should be far enough from the maximum of the distribution that neglecting the 
(small) number of particles in the boundary bins does not introduce a noticeable error. Requiring 
the total number of electrons to not change with time, one obtains certain conditions on the 
differencing scheme at the neighborhood of the boundaries. These conditions lead us to choose the 
differencing scheme for the second point on the energy axis to be 



fn+l 
J2 
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and for the next to last point, 

rn+l _ rn I / "m/m + "M-i/m-A . \ r,n fn , j^n fn ] f AK\ 

}m - /m + At AxM-i + AxM / 2 ^ [-^m}m + ^m-i/m-iJ , (^5) 

where M = N — 1, and A'' is the total number of points in energy space. For the third point (0:3), 
the only change needed to be made is to set the coefficient p in Equation (Al) to p = q. 

Because Equations (A4) and (A5) employ only two boundary points instead of three as any 
other A; 7^ 2, M , it is not possible to make the scheme conserve particle number and energy 
simultaneously for these two bins. We need two parameters to conserve these integrals, while we 
have only one parameter available in Equations (A4) and (A5) (i.e., the ratio a/ p). Physically, by 
bouncing particles scattering into the first bin back to the second and third bins {fi{t) = at all 
times) , we conserve the number of particles, but change the energy of these particles to make them 
stay inside this energy range (i.e., inside the region X2 < 7 — 1 < Xjv-i)- This process therefore 
does not conserve energy, and we must subtract this excess energy from the total energy change 
during each iteration. Practically, this is done by introducing a small correction to the electron 
energy exchange coefficient, which allows us to conserve energy up to double precision accuracy 
during the entire particle thermalization time. 

Note that the mirror boundary conditions are perfectly physical and are applicable to any 
electron distribution function with large or small deviations from a Maxwellian. Indeed, if there 
were a physically substantial flow through the low energy boundary, that would mean that there 
is either a source or sink of particles with E <^ Q^,, where Qg is the dimensionlcss electron 
temperature, Qg = kTe/nieC^. While this situation is mathematically possible, it is not plausible 
physically (if there is an influx of cold electrons from outside into the system, one can always choose 
the low energy boundary such that it includes these electrons). On the high energy end, any electron 
acceleration mechanism can only accelerate electrons to some maximum energy (specific for this 
mechanism), after which electron cooling overcomes the heating. Therefore, if one chooses the high 
energy boundary much higher than this maximum energy, then there can be no flux through the 
boundary, since the electrons get 'turned back' to lower energies before they reach the boundary. 

The ultimate test of the given scheme is to see how close the equilibrium distribution function 
approaches a Maxwellian (when one includes only the lepton-lepton FP part of the full lepton 
Boltzmann equation). We have run a number of tests (described in the next few sections) and 
find excellent agreement between the equilibrium distribution and a Maxwellian (with a relative 
deviation of less than 2%). 

A. 2 Numerical Stability Of The Scheme 

Equation (Al) is a non-linear integro-differential equation, for which strictly speaking, the 
usual Von Neumann analysis should not be used to test its stability. Fortunately, the fact 
that the FP coefficients are integrals over the distribution function makes the stability analysis 
essentially linear. Conducting the stability analysis of our scheme as if a{E,t) and D{E,t) are 
fixed functions, or at most, slowly varying functions compared to f{E,t), and using the usual 
Von-Neuman approach (e.g.. Press et al. 1986), we arrive at the approximate stability criterion: 

r/ D^At \2 akAt 2\ ^ . . ... 
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Equation (A6) is a necessary condition for the scheme (Eq. Al) to be numerically stable. While 
this condition is only approximate, we find that by choosing a time step roughly equal to or smaller 
than that which satisfies Equation (A6), our differencing scheme is always stable. 

For applications where 'exact' energy conservation is not an issue (e.g., in systems with a low 
optical depth and a high compactness, such that the energy is radiated efficiently from the system), 
we were able to develop a half implicit-half explicit scheme that is much faster than the explicit 
scheme (see the next section). 

A.3 The Numerical Procedure 

We follow the procedure outlined in §A.l above, except that we here include the additional 
terms corresponding to Compton interactions and proton heating. The particle distributions are 
evolved more often than the photon distribution, i.e., the photon distribution is kept fixed during 
a time AtM = MAt while the particles proceed through their evolution a number M (s> 1) of At 
steps. This approach is necessary because otherwise in each time advancement, one needs to re- 
evaluate the photon annihilation cross section and the Compton scattering matrix, among others. 
This increases the computation time considerably. The value of At is dictated by the stability 
criterion given in Equation (A6). It turns out that this condition is much stricter than any other 
criterion for At based, e.g., on the requirement that distribution functions change very little during 
one time advancement (i.e. |A/(£^)| ^ \f{E)\). The latter should be used for determining the 
time step when evolving the photon distribution since there is no numerical instability problem 
analogous to Equation (46). We are therefore justified in using AtM 3> At. 

The computation time depends considerably on whether one uses explicit or implicit schemes 
and on the Thomson optical depth of the system. For a Thomson optical depth rr ~ 1, the explicit 
scheme takes a few days on a Pentium PC running Linux (for 70 bins in the electron energy space 
and 70 bins for photons) , whereas an implicit scheme takes from 6 hours to a day. Obviously, such 
a long computing time makes it difficult to use the described methods in practice, for example to 
fit observations. 

For most of the applications, however, a time-independent code should suffice. We have been 
able to write and use such a code successfully. Using the approximate solution of the full electron 
equation given by Equation (50), and assuming that the positron distribution is same as that of 
the electrons, we were able to reduce computing time to 2 — 5 minutes. In essence, the code is 
not any slower than a corresponding code that assumes that particles are thermal and takes into 
account the same interactions and follows radiation transport in the same manner. Thus, the 
time-independent code is highly efficient and might easily be used in X-ray fitting packages, such 
as XSPEC. This code and some results are to be published elsewhere (the results presented in the 
current paper were obtained using the full time-dependent implicit code) . The authors can provide 
their code for use by anybody who feels a need to use exact particle distributions. 

To veriiy that our code was correctly resolving the particle distribution in energy space, we also 
varied the number of energy points. A comparison of our results from a test with 100 energy bins 
(for both electrons and photons) with those from a simulation with 70, showed that the relative 
differences in the distribution, spectrum, temperature and the number of positrons were ^ 1%. 
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A.4 Large Angle Scatterings 

As already mentioned in §3, it has been known since the early work of Chandrasekhar and 
Landau that large angle scattering events are not important for Coulomb interactions in a gas, since 
the rate of energy exchange due to these is much smaller than that due to small angle collisions. 
This point is worth discussing further here, since it appears to be the basis for the difference 
between our approach and that of Dermer & Liang 1989. Large angle scatterings correspond to 
small impact parameters. We can estimate the latter for large angle scatterings r; by writing 

— > lAElmeC^ (A7) 

n 

(recalling that = 7 — 1 is dimensionless) . The cross section is then ai ( A£^) ^ irrf , and the energy 
exchange rate is 

ai{E) ~ ai{AE)\AE\n,pc < tt n,c(3/\AE\ . (A8) 

Let us consider the high energy portion of the electron distribution, i.e., E Q. Large angle 
scatterings correspond to a large change in the electron energy AE ~ E (Equation [14]). Therefore, 
we can put AE = E in the above equation, and re- write it as 

Equation (30) can be used to estimate the small angle energy exchange rate 

where (7) is the average thermal 7-factor. The ratio of the large angle energy exchange rate to 
that of small angle scatterings is thus 

^4^.^ 1 «i, (All) 

as{E) 7-1 41nA ' ^ ^ 

as long as InA S> 1. Let us now assume E <^ Q. In this case we should use the average thermal 
P and 7 in Equation (A9), since the energy exchange can be as large as ~ 0, and /3 in equation 
(AS) is really the relative velocity of the two colliding particles, which is of the same order as the 
thermal velocity. The ratio of the scattering rates becomes 

ai(E) (/3)27 11^ ^ ^ 

^ ~ ' < . < 1 (A12) 

a,(E) (7)-l 41nA 41nA ^ ^ 

A similar estimate holds for energies E ~ (E), except for the narrow region where the small angle 
energy-exchange coefficient goes to zero (Fig. [3]). 

Therefore, it is evident that one may neglect the large angle scatterings in the full kinetic 
equation that describes Coulomb collisions. Moreover, we argue that it is even necessary to do so 
in the calculation of the FP coefficients, since otherwise the approximation made for small scattering 
angles becomes invalid for large angles and can lead to wrong conclusions. Physically, one cannot 
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decompose the electron distribution function in the original Boltzmann kinetic equation, since the 
large angle scatterings do not constitute a diffusion process. For example, the use of the diffusion 
coefHcient as computed by Dermer & Liang (1989) can lead to a power-law electron distribution 
on the high energy end instead of a Maxwellian distribution, if other interactions, such as inverse 
Compton cooling are not important. The large angle scattering terms, which according to the 
above discussion should contribute very little, in fact start to dominate the diffusion coefficient. 
The energy exchange coefficient, on the other hand, is unaffected by the inclusion of these terms 
since (by coincidence) averaging over 4>* cancels out the large positive and negative AE in Equation 
(14). 
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